﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace NCS_CS
{
    class Lfit
    {
        const double TINY = 1.0e-27;


        //static public double poly(double x, int grado, double[] coeff)
        //{
        //    double p = coeff[grado];
        //    for (int i = grado - 1; i >= 0; i--)
        //    {
        //        p = p * x + coeff[i];
        //    }
        //    return p;
        //}



        static public double[][] matrix(int nrow, int ncol)
        {
            double[][] m;

            m = new double[nrow][];
            for (int i = 0; i < nrow; i++)
                m[i] = vector(ncol);

            for (int i = 0; i < nrow; i++)
                for (int e = 0; e < ncol; e++)
                    m[i][e] = 0;

            return m;
        }

        static public double[] vector(int nrow)
        {
            double[] v;

            v = new double[nrow];

            return v;
        }

        //double[][] matrixl(int nrow, int ncol)
        //{
        //    double[][]m;

        //    m = new double[nrow][];
        //    for (int i=0; i<nrow; i++)
        //        m[i] = vectorl(ncol);

        //    for (int i=0; i<nrow; i++)
        //        for (int e=0; e<ncol; e++)
        //            m[i][e] = 0;

        //    return m;
        //}

        //long double* vectorl(int nrow)
        //{
        //    long double* v;

        //    v = new long double[nrow];

        //    return v;
        //}
        //void test(int n)
        //{
        //}

        static public void lfit(double[] x, double[] y, int[] p, int ndat, ref double[] a,
            int ma, ref double chisq)//, void (*funcs)(double, double [], int))
        {
            int i, j, k, l, m;
            double ym, wt, sum;
            double[][] beta;
            double[] afunc;
            double[][] covar;
            if (ma == 0)
                return;
            beta = matrix(ma, 1);
            afunc = vector(ma);
            covar = matrix(ma, ma);

           
            //::MessageBox.Show("LFIT: no parameters to be fitted");

            for (j = 0; j < ma; j++)
            {
                for (k = 0; k < ma; k++)
                    covar[j][k] = 0.0;
                beta[j][0] = .0;
            }
            for (i = 0; i < ndat; i++)
            {
                dpoly(x[i], afunc, ma);
                ym = y[i];

                for (j = -1, l = 0; l < ma; l++)
                {
                    wt = afunc[l] * p[i];
                    for (j++, k = -1, m = 0; m <= l; m++)
                        covar[j][++k] += wt * afunc[m];
                    beta[j][0] += ym * wt;
                }
            }

            for (j = 1; j < ma; j++)
                for (k = 0; k <= j; k++)
                    covar[k][j] = covar[j][k];
            gaussj(covar, ma, beta, 1);
            for (j = -1, l = 0; l < ma; l++)
                a[l] = beta[++j][0];
            ///////////////result//ma

            chisq = 0.0;
            for (i = 0; i < ndat; i++)
            {
                dpoly(x[i], afunc, ma);
                for (sum = .0, j = 0; j < ma; j++)
                    sum += a[j] * afunc[j];
                chisq += System.Math.Pow(y[i] - sum, 2.0);
            }

        }


        static public void gaussj(double[][] a, int n, double[][] b, int m)
        {
            int[] indxc, indxr, ipiv;
            int i, icol = 0, irow = 0, j, k, l, ll;
            double big, dum, pivinv;//, temp;

            indxc = new int[n];
            indxr = new int[n];
            ipiv = new int[n];

            for (j = 0; j < n; j++)
                ipiv[j] = -1;

            for (i = 0; i < n; i++)
            {
                big = .0;
                for (j = 0; j < n; j++)
                    if (ipiv[j] != 0)
                        for (k = 0; k < n; k++)
                        {
                            if (ipiv[k] == -1)
                            {
                                if (System.Math.Abs(a[j][k]) >= big)
                                {
                                    big = System.Math.Abs(a[j][k]);
                                    irow = j;
                                    icol = k;
                                }
                            }
                            else
                                if (ipiv[k] > 0)
                                    return;
                        }
                ++(ipiv[icol]);
                if (irow != icol)
                {
                    for (l = 0; l < n; l++)

                        SWAP(ref a[irow][l], ref a[icol][l]);
                    for (l = 0; l < m; l++)
                        SWAP(ref b[irow][l], ref b[icol][l]);
                }
                indxr[i] = irow;
                indxc[i] = icol;
                if (a[icol][icol] == .0)
                    return;
                pivinv = 1.0 / a[icol][icol];
                a[icol][icol] = 1.0;
                for (l = 0; l < n; l++)
                    a[icol][l] *= pivinv;
                for (l = 0; l < m; l++)
                    b[icol][l] *= pivinv;
                for (ll = 0; ll < n; ll++)
                    if (ll != icol)
                    {
                        dum = a[ll][icol];
                        a[ll][icol] = .0;
                        for (l = 0; l < n; l++)
                            a[ll][l] -= a[icol][l] * dum;
                        for (l = 0; l < m; l++)
                            b[ll][l] -= b[icol][l] * dum;
                    }
            }

            for (l = n - 1; l >= 0; l--)
            {
                if (indxr[l] != indxc[l])
                    for (k = 0; k < n; k++)
                        SWAP(ref a[k][indxr[l]], ref a[k][indxc[l]]);
            }

        }

        static void dpoly(double x, double[] p, int np)
        {
            p[0] = 1.0;
            for (int j = 1; j < np; j++)
                p[j] = p[j - 1] * x;
        }


        //static public void mrqmin(double[] x, double[] y, double[] sig, int ndata, double[] a,
        //            int ma, double[][] covar, double[][] alpha, ref double chisq,
        //            myfuncs funcs, ref double alambda)
        //{
        //    int j, k, l;
        //    int mfit = 0;
        //    double ochisq = 0;
        //    double[] atry = null, beta = null, da = null;
        //    double[][] oneda = null;

        //    if (alambda < .0)
        //    {
        //        atry = vector(ma);
        //        beta = vector(ma);
        //        da = vector(ma);
        //        mfit = ma;
        //        oneda = matrix(ma, 1);
        //        alambda = .001;
        //        mrqcof(x, y, sig, ndata, a, ma, alpha, beta, ref chisq, funcs);
        //        for (j = 0; j < ma; j++)
        //            atry[j] = a[j];
        //        ochisq = chisq;
        //    }

        //    for (j = 0; j < mfit; j++)
        //    {
        //        for (k = 0; k < mfit; k++)
        //            covar[j][k] = alpha[j][k];
        //        covar[j][j] = alpha[j][j] * (1.0 + alambda);
        //        oneda[j][0] = beta[j];
        //    }

        //    gaussj(covar, mfit, oneda, 1);

        //    for (j = 0; j < mfit; j++)
        //        da[j] = oneda[j][0];

        //    //if (*alambda == .0)
        //    //{
        //    //    free_matrix(oneda,ma);
        //    //    delete [] da;
        //    //    delete [] beta;
        //    //    delete [] atry;
        //    //    return;
        //    //}

        //    for (l = 0; l < ma; l++)
        //        atry[l] = a[l] + da[l];

        //    mrqcof(x, y, sig, ndata, atry, ma, covar, da, ref  chisq, funcs);

        //    if (chisq < ochisq)
        //    {
        //        alambda *= .1;
        //        ochisq = chisq;
        //        for (j = 0; j < mfit; j++)
        //        {
        //            for (k = 0; k < mfit; k++)
        //                alpha[j][k] = covar[j][k];
        //            beta[j] = da[j];
        //        }
        //        for (l = 0; l < ma; l++)
        //            a[l] = atry[l];
        //    }
        //    else
        //    {
        //        alambda *= 10.0;
        //        chisq = ochisq;
        //    }

        //}
        //public delegate void myAction(List<Cmd> i);
        //public delegate void myfuncs(double a, double[] b, ref double c, double[] d, int e);

        //static public void mrqcof(double[] x, double[] y, double[] sig, int ndata, double[] a,
        //            int ma, double[][] alpha, double[] beta, ref double chisq,
        //            myfuncs funcs)
        //{
        //    int i, j, k, l, m;
        //    double ymod = 0, wt, dy;
        //    double[] dyda;
        //    double sig2i;

        //    dyda = vector(ma);

        //    for (j = 0; j < ma; j++)
        //    {
        //        for (k = 0; k <= j; k++)
        //            alpha[j][k] = .0;
        //        beta[j] = .0;
        //    }

        //    chisq = .0;

        //    for (i = 0; i < ndata; i++)
        //    {
        //        funcs(x[i], a, ref ymod, dyda, ma);
        //        sig2i = 1.0;///(sig[i]*sig[i]);
        //        dy = y[i] - ymod;
        //        for (l = 0; l < ma; l++)
        //        {
        //            wt = dyda[l] * sig2i;
        //            for (m = 0; m <= l; m++)
        //                alpha[l][m] += wt * dyda[m];
        //            beta[j] += dy * wt;
        //        }
        //        chisq += dy * dy * sig2i;
        //    }

        //    for (j = 1; j < ma; j++)
        //        for (k = 0; k < j; k++)
        //            alpha[k][j] = alpha[j][k];

        //    //delete []dyda;
        //}

        //static public void prova(double x, double[] a, ref double y, double[] dyda, int na)
        //{
        //    int i;
        //    y = 0.0;
        //    for (i = na - 1; i >= 0; i--)
        //    {
        //        y = y * x + a[i];
        //    }

        //    dyda[0] = 1;
        //    for (i = 1; i < na; i++)
        //        dyda[i] = dyda[i - 1] * x;
        //}

        //static public void retta(double x0, double y0, double x1, double y1, ref double a, ref double b)
        //{
        //    b = (y1 - y0) / (x1 - x0);
        //    a = y0 - (b) * x0;
        //}


        //static public void pearsn(double[] x, double[] y, int n, ref double r)
        //{
        //    int j;
        //    double yt, xt;
        //    double syy = 0.0, sxy = 0.0, sxx = 0.0, ay = 0.0, ax = 0.0;

        //    for (j = 0; j < n; j++)
        //    {
        //        ax += x[j];
        //        ay += y[j];
        //    }
        //    ax /= n;
        //    ay /= n;
        //    for (j = 0; j < n; j++)
        //    {
        //        xt = x[j] - ax;
        //        yt = y[j] - ay;
        //        sxx += xt * xt;
        //        syy += yt * yt;
        //        sxy += xt * yt;
        //    }
        //    r = sxy / (System.Math.Sqrt(sxx * syy) + TINY);
        //}

        //static public void avevar(double[] data, int n, ref double ave, ref double var)
        //{
        //    int j;
        //    double s, ep;

        //    for (ave = 0.0, j = 0; j < n; j++)
        //        ave += data[j];
        //    ave /= n;
        //    var = ep = 0.0;
        //    for (j = 0; j < n; j++)
        //    {
        //        s = data[j] - (ave);
        //        ep += s;
        //        var += s * s;
        //    }
        //    var = (var - ep * ep / n) / (n - 1);
        //}


        //static double round(double c)
        //{
        //    //double result;
        //    //double i;
        //    //double d = modf(c, &i);
        //    //if (d > .5)
        //    //    result = ceil(c);
        //    //else
        //    //    result = floor(c);
        //    return c;
        //}
        static public void SWAP(ref double a, ref double b)
        {
            double tmp;
            tmp = a;
            a = b;
            b = tmp;
        }
    }
}
